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ABSTRACT 


Flow fields involving turbulent flames in premixed gases under a 
variety of conditions are modeled by the use of a numerical technique 
based on the random vortex method to solve the Navier-Stokes equations and 
a flame propagation algorithm to trace the motion of the front and imple- 
ment the Huygens principle, both due to Chorin. A successive over-relaxa- 
tion hybrid method is applied to solve the Euler equation for flows in an 
arbitrarily shaped domain. The method of images, conformal transformation, 
and the integral-equation technique are also used to treat flows in special 
cases, according to their particular requirements. 

Salient features of turbulent flame propagation in premixed gases are 
interpreted by relating them to the aerodynamic properties of the flow 
field. Included among them is the well-known cellular structure of flames 
stabilized by bluff bodies, as well as the formation of the characteristic 
tulip shape of flames propagating in ducts. In its rudimentary form, the 
mechanism of propagation of a turbulent flame is shown to consist of (i) 
rotary motion of eddies at the flame front, (ii) self-advancement of the 
front at an appropriate normal burning speed, and (iii) dynamic effects of 
expansion due to exothermicity of the combustion reaction. An idealized 
model is used to illustrate these fundamental mechanisms and to investigate 
basic aerodynamic features of flames in premixed gases. 

The case of a confined flame stabilized behind a rearward-facing step 
is given particular care and attention. Solutions are shown to be in satis- 
factory agreement with experimental results, especially with respect to 
global properties such as the average velocity profiles and reattachment 
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length. Velocity fluctuations are found to compare quite well with experi- 
mental data, exhibiting discrepancies which can be plausibly ascribed to 
three-dimensional effects and the scarcity of the numerical data sample. 
The unconfined turbulent flow behind a circular flameholder is calculated 
to study the detailed mechanism of flame stabilization. The process of 
vortex shedding and the development of the Von Karman vortex street at a 
high Reynolds number are described in terms of the vorticity field and the 
instantaneous streamline pattern. Although the flame front is stabilized on 
the outer contours of eddies, the classical alternating large scale struc- 
ture of the wake is destabilized by the expansion due to the exothermicity 
of combustion. The parametric study of flame propagation indicates that the 
flame spread and the burning rate are relatively insensitive to the 
Reynolds number in the range under consideration, whereas the product of 
the normal burning speed and the normalized density ratio increment across 
the front, expressing the action of the flame as a velocity source, is of 
crucial importance. 
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Nomenclature 


: Area 

a : coordinate of a point on the real axis of the f - plane 
a : radius of a cylinder 

b : coordinate of a point on the real axis of the f - plane 

c : coordinate of a point on the real axis of the f - plane 

D : molecular diffusivity 

d : diameter of a cylinder 

F : differential transformation function 

f : volume fraction of burned gases in a cell 
G : Gaussian distribution function 
H : width of the duct 

h : mesh size 

I : number of grid points in x - direction 
i : V— 1 

J : number of grid points in y - direction 
K : coefficient of transformation function 

k : thermal conductivity 

k : turbulence kinetic energy 

Le : Lewis number, defined as D / k 
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M : total number of grid points for finite-difference calculation 
N : total number of vortex elements 
n : unit normal vector 

P : pressure 

p : coordinate of a point on the surface of a bluff-body 

q : coordinate of a point outside the body 

Re : Reynolds number, defined as Ul/v 

rj : position vector of flame front 

r 0 : core radius of a vortex blob 

r, : core radius of a source blob 

Su : laminar flame speed 

s : unit vector tangential to streamline 

T : nondimensional time 

A t : nondimensional time step 

U : mean velocity 

u : velocity vector 

u : velocity component in x - direction 

V : density ratio across flame front 

v : velocity component in y - direction 

W : complex velocity 

x : x - coordinate in the physical plane 

y : y - coordinate in the physical plane 
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z : = x + i y , complex coordinate 

a : an interior angle of a polygon in the physical plane 

fl : an interior angle of a polygon in the physical plane 

T : strength of a vortex 

7 : an interior angle of a polygon in the physical plane 

A : strength of a source blob 

5 : boundary-layer thickness 

e : local rate of expansion 

£ : complex coordinate in the transformed plane 

T) : random displacement 

: angle measured counter-clockwise from x - axis to the tangent of the 
wall 

X : a parameter used to determine the over-relaxation coefficient 
v : kinematic viscosity 
£ : vorticity 

p : density 

a : standard deviation 

$ : velocity potential 

: equivalence ratio 
't' : stream function 

o : over-relaxation coefficient 

subscript 
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b : burned medium 

c : combustion 

f : flow field 

i : x index of grid point 

j : y index of grid point 

k : index of grid point along the wall 
max : maximum value 
n : normal component 

p : potential component 

s : effect due to source blobs 

u : unburned medium 

v : effect due to vortex blobs 

w : value at wall 
eo : value at infinity 

superscript 

n : index of time step 
• : conjugate of a complex number 

o : degree 
’ : fluctuation 
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CHAPTER 1 


Introduction 

Turbulent reacting flows have attracted lots of attentions in recent 
years. This is because that they not only involve all kinds of chemical reac- 
tions, energy transports and complicated fluid mechanical turbulence, thus 
making them challenging, but also because there are still a number of 
unresolved problems in this field due to the complicities. In premixed tur- 
bulent combustion, the thermal conduction and diffusion from the hot pro- 
ducts to preheat the reactants can be augmented both indirectly by distor- 
tion of flame surfaces and directly by turbulent mixing. The mean reaction 
rate is thus more strongly influenced by the turbulence than by chemical 
kinetics factor, so that premixed turbulent combustion is primarily a com- 
plex fluid mechanical problem. 

Ever since Damkohler's (1940) classical theoretical and experimental 
study of turbulent flames in premixed gases, the wrinkle laminar flame idea, 
which considers large-scale turbulence as simply distorting the flame 
without affecting its structure or propagating speed, has been playing a very 
important part in today's research of turbulent flame. When the radius of 
curvature of a steady distorted flame is considerably larger than the width 
of the thermal structure of the flame, a flame model which simply treated 
the flame as a hydrodynamic discontinuity of density and assumed that the 
propagating speed relative to the gas along the flame did not vary, has been 
proposed (Landau, 1944, Sivashinsky, 1978), the latter studies a number of 
flame problems with great deal of success. 
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Vortex method was used to calculate inviscid flows at the early stage 
(Rosenhead, 1931, Takagi, 1984, and Clements, 1973). Later, it was 
developed (Chorin, 1973) to include the viscous effect by vorticity genera- 
tions at boundaries and random walks to simulate the diffusion process in a 
statistical sense. Without resorting to closure modeling of turbulence and 
introducing numerical diflusivity by traditional finite-difference method, the 
random vortex method is able to compute the evolution and interaction of 
he turbulent flow field at high Reynolds number. It has thus been used to 
calculate the two-dimensional turbulent flows quite successfully (Ashurst, 
1979, Cheer, 1982, Ghoniem et al, 1982, Chen et al, 1983 and Hsiao et al 
1984). 

Since the diffusion equation is adroitly handled by random walks of vor- 
tex elements, the remained problem is the solving of the Euler's equation. 
The vorticity of the flow field is represented by a number of vortex blobs. 
There are many ways of finding the velocity field induced by vortex blobs in 
a specific domain. Some examples of the methods used to handle the zero 
normal velocity condition will be outlined and discussed in our numerical 
model. Their applicability to different flow configurations and evaluations of 
each method will also be given. 

To preserve the grid-free principle of the vortex method, the method of 
images, conformal transformation and the integral-equation method can be 
applied to solve the inviscid flow field without imposing grids in the computa- 
tional domain. The velocity fields are found by direct interactions of all vor- 
tex blobs. The computational time is thus proportional to N , where N is the 
number of vortex elements. For most practical cases, transformation func- 
tions have to be evaluated by complicated numerical integration and the 
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source-density functions of the integral-equation method along body sur- 
faces also have to be determined by solving a set of linear algebraic equa- 
tions. These not only bring additional difficulties and errors into the numer- 
ical method, but also increase computational efforts dramatically. 

A most straightforward way to deal with all kinds of geometries is the 
finite-difference method. The so-called fast Poisson solvers (Buzbee et al, 
1970, Swarztrauber et al, 1975) are frequently used to find the stream func- 
tions or velocity potentials at grid points. The velocities can then be 
obtained by differentiation and interpolation. Since the existing numerical 
codes in solving Poisson equation can not be very generally applied to any 
flow domain, the successive over- relaxation method is suggested here to 
solve the Poisson equation in arbitrary- shape domain because of its easiness 
of programming and the competitive efficiency. 

Recently, Oppenheim and Ghoniem (1983) were able to deduce three 
fundamental mechanisms of turbulent flame propagation, including rotary 
motion of eddies with different length scales, flame advancement in the nor- 
mal direction, and volumetric expansion due to exotbermicity, from many 
interesting flame phenomena in premixed combustion. These three mechan- 
isms are also the basic components governing flame propagation in our 
numerical model. The interaction between flame fronts and turbulent 
eddies can result in the distortion of the flame, so that the surface area is 
markedly increased. On the other hand, the entrainment mechanism of the 
recirculation zone or large-scale structures can enhance the mixing by 
bringing more fresh reactive mixtures into products. The result is an 
increase in the reaction rate. While the volumetric expansion due to 
combustion reaction not only modifies the flow field, such as the reduction 



of the recirculation zone and the divergence of streamlines outside flame 
boundaries, but also changes the flame shape by the feedback action of 
local velocity. 

Flame propagation in a closed volume has somehow different natures 
from that in an open volume. The former is always accomplished by a con- 
tinuous rise of system pressure. This will cause a temperature rise in the 
mixture so that the reaction rate or flame propagation speed also increases. 
The compression of the reactive mixtures and the flow field generated by 
the combustion can result in several very interesting phenomena of flame 
deformation. The cause of the so-called "tulip” shape flame in an enclosed 
duct is still an unresolved problem. By assuming a slowly propagating flame 
so that the pressure rise in the enclosed vessels can be instantaneously bal- 
anced over the entire volume, Sivashinsky (1979) was able to derive a hydro- 
dynamic theory to describe a centrally-symmetric flame in an enclosed 
volume. Majda and Sethian (1984) also accomplished a theory of zero- 
Mach-number combustion using a different approach. They applied it to 
study the flame propagation of the swirling flow inside a square vessel. 

Unfortunately, these theories have not been able to apply to problems 
which will demonstrate certain interesting features of flame propagation in 
a closed domain, such as the deceleration of the flame motion, the deforma- 
tion of flame shape and the detailed flow field result from the reaction. 
Before proceeding to a more complete theory or a detailed numerical model, 
we will use an ideal potential flow model in Chapter 3 to present , we think, 
the essence of the cause of the tulip-shape flame from a purely hydro- 
dynamic point of view. The individual or combined effects of the three fun- 
damental mechanisms of flame propagation will be used to illustrate the 



interpretation. 

The turbulent shear layer behind a rearward-facing step in a channel 
has been extensively investigated for both reacting and non-reacting flows 
in connection with their significance for a number of practical engineering 
systems and the simplicity of the case where the separated flow region is 
fixed. Summaries of experimental studies of non-reacting flows with this 
specific geometry can be found m the review papers by Eaton & Johnson 
(1981) and Durst 8c Tropea (1982). Some numerical solutions of this flow 
configuration can be also found in the publications of Durst & Rastogi 
(1979), Ashurst (1979), Ghoniem et al (1982), Dai et al (1983) and Walterick 
et al (1984). Many things have been learned about the two-dimensional 
mechanism of flow separation and reattachment. The existing data of mean 
flow properties and turbulence characteristics of such flows, although still 
spread out a finite range, are sufficiently satisfactory to compare with other 
kinds of flows and to be used for checking numerical models and computer 
codes. 

The study of the reacting flow for a rearward-facing step combustor was 
first performed by Ganji 8c Sawyer (198D) using high speed schlieren photog- 
raphy. They discovered that the reacting flow is dominated by large scale 
coherent structures in mixing layers. Pitz & Daily (1983) then measured the 
detailed flow properties and turbulence statistics up to the fourth moment. 
They found a shift of the peak positions of the turbulence intensity profiles 
toward the recirculation zone and a reduction of the reattachment length 
for reacting flow. The streamwise velocity PDFs were also shown to be 
double-peaked in the reacting shear layer near the step. Meanwhile, 
Ghoniem et al (1982) developed a numerical model using random vortex 
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method and a flame propagation algorithm by Chorin (1980) to simulate the 
turbulent flow field behind the rearward-facing step. They obtained a quali- 
tative description of the large scale structure of the turbulent flow as well 
as the salient feature of the flame front in agreement with the experimental 
results obtained by Ganji & Sawyer. 

Since the upstream flow conditions have been found to have a 
significant effects on the flow development behind the step, the computa- 
tional domain is extended to cover the contraction inlet section for a more 
complete simulation of the experiment. The mean flow velocities and tur- 
bulence intensities will be computed in addition to the development of vorti- 
city field and flame fronts. The comparison between numerical results and 
experimental data will also be carried out. 

A circular cylinder is the simplest geometry of bluff bodies. The investi- 
gation of the wake behind a cylinder have been done on flows at different 
Reynolds numbers (Taneda, 1952, Roshko, 1981, Gerrard, 1988, Mair & Maull, 
1971 and Perry et al, 1982). A recent experimental study by Cantwell & Cole 
(1983) gave detailed entrainment and transport processes in the near wake 
region at very high Reynolds number. However, the vortex-shedding 
mechanism in the formation zone behind the cylinder is still not well under- 
stood. 

Flame stabilized by a rod holder is also an extensively studied subject 
(Peterson & Emmons, 1981, Cheng 8t Ng, 1983, Katsuki et al, 1983 and 
Yoshida & Tsuji, 1984) in combustion area. Beside all the turbulence statis- 
tics, the flame stability and turbulence-combustion interaction seem to be 
the two most frequently studied topics. It was conjectured that the vortex- 
shedding process in the near wake region triggers the wrinkles of flame 
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front downstream. 


The connection between the study of non-reacting flow and reacting 
flow in the wake region is poor. Owing to the lack of proper flame model, 
most research works of this subject have been concentrating on the experi- 
mental measurement of velocity, temperature and density fluctuations (Kil- 
ham Sc Kirmani, 1979, Susuki et al, 1979 and Dandekar et al, 1982). Their 
results, although quantitative in nature, are insufficient to the understand- 
ing of the mechanism of flame wrinkling and turbulence-combustion 
interaction. It has been suggested that the qualitative features of the flame 
shown by flow visualization would be useful to explore the possible interrela- 
tionship between the flame structures and the observed behavior of the tur- 
bulence statistics and correlations. 

Our numerical model, being capable of simulating large scale structures 
in the flow field and providing a phenomenological description of the flame 
shape and flame-flow interaction, seems to be a proper tool in serving a link 
between the quantitative measurements and the observed flame structures. 
In Chapter 5 the detailed investigation of the flow development and vortex- 
shedding process for non-reacting flow is conducted, first, to demonstrate 
the flow phenomena in the wake region. The flame front which starts from 
the rear stagnation point of the cylinder follows the motion of turbulent 
eddies and advances into its normal direction in the meantime. The 
volumetric expansion, which is a major factor of modifying the flow field, is 
also counted. The resulting flame front and flow field due to the combined 
effect of these three fundamental mechanisms can certainly let us proceed 
to the further understanding of the unconflned flame stabilized by a circu- 
lar cylinder. 
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Chapter 2 


Numerical Uodel 


The numerical model used here is in principle based on that developed 
by Ghoniem et al (1982) in the paper studying turbulent flow in a combus- 
tion tunnel. The model uses random vortex method to calculate turbulent 
flows and is currently applied to two-dimensional cases only, thus an impor- 
tant vorticity-maintenance mechanism in turbulence known as vortex 
stretching is absent in our analysis. To investigate the aerodynamic features 
of the turbulent combustion, the flame front is treated as an infinitely thin 
interface separating incompressible premixed reactants from the products. 
The normal burning speed along flame front is assumed to be as a 
prescribed constant. Effects of chemical kinetics and flame structure are 
therefore not taken into account. Turbulent transport and flame propaga- 
tion are linked together by the velocity field. Flame front is convected by 
the local velocity while the dynamic effect due to exothermicity of combus- 
tion is modeled by volumetric sources located along flame fronts according 
to amount of burning to augment the flow field. 


2.1. Turbulent Non-reacting Flows 

Turbulent flows have been studied for more than a century, but no gen- 
eral approach to the solution of problems in turbulence exists (Tennekes & 
Lumley, 1970). However, equations which are used to describe the tur- 
bulence in fluids are generally accepted to be continuity equation and 
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Navier-Stokes equation. In two-dimensional incompressible fluids they can 
be expressed in the non-dimensional form a9 

V ■ u = 0 (2.1) 

= -7p + Re -1 7 8 u (2.2) 

with boundary condition 

u = 0 (2.3) 

where u is velocity, p is pressure and Re is Reynolds number, while D/Dt is 
material derivative, 7 is del operator and 7^ is the Laplace operator. 

In most turbulent-flow calculations, the Reynolds number is usually so 
large that traditional flnite-diflerence method is difficult to apply especially 
near boundaries where velocity gradients are large. The random vortex 
method which involves the simulation of the process of vorticity generation 
and disposal, developed by Chorin (1973), has been found to be able to over- 
come these problems and calculates turbulent flows of high Reynolds 
numbers with plausible accuracy without resorting to closure models. 
Governing equations under consideration are vorticity transport equations 

|f = Re' 1 V 8 ( (2.4) 

and 

T 8 * = - ( (2.5) 


3* 3* 

u = — ; v - 

dy dx 


where £ = 7 x u is vorticity and ^ is stream function. 


( 2 . 6 ) 


The solution procedure of above equations is based on fractional steps, 
according to which the vorticity equation can be split into a sum of convec- 
tion component (the Euler equation) and diffusion component (diffusion 
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equation). The result obtained by adding up solutions to these two equa- 
tions in succession has been proven to converge by Chorin et al (1978) and 
Held (1979). 


2.1.1. Solutions of the Euler Equation 

Inviscid flow field is solved by discretizing the vorticity field into a 
number of vortex blobs. To get rid of the singularity at the center of a point 
vortex, the vorticity of each blob is constructed to have small but finite sup- 
port. The motions of vortex blobs are described by the inviscid Euler equa- 
tion 

0 or + (u ■ V ) f = 0 (2.7) 

and 

V z * = - £ (2.8) 

In an infinity flow domain without any boundary, the velocity field can 
be easily found by summing up all velocities induced by vortex blobs accord- 
ing to Biot-Savart's law 


Uv - i v v 



'?J ..... 

Z - Zj 

1 

max(| 

Z — Zj 

. rj ( * “ ) 


(2.9) 


where i = (-1)^ , z = x + iy and r Q , Tj are the core radius and the strength 
of the vortex blob respectively. In most practical cases, however, rigid walls 
exist in the flow domain. The zero normal velocity condition 


u ■ n = 0 

has to be satisfied at boundaries. There are several different ways of solving 
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the above equations. Four examples of these methods which will be used 
together with vortex method to solve the inviscid flow are presented in fol- 
lowing sections. 

2.1. 1.1. Method of Images 

If the geometry of the boundary is simple such as an infinity plate or a 
circular cylinder, the inviscid flow solution can be obtained by the method of 
images. The solution of a uniform potential flow past a circular cylinder can 
be found by the circle theorem (Milne-Thomson, 1971). For every vortex 
blob in the flow domain, there will be a corresponding vortex blob with equal 
strength and opposite sign at its image location. The normal velocity com- 
ponent induced by each pair of vortices at the boundary can thus cancel 
out each other. Inviscid flow solutions obtained by the use of this method 
always turn out to be exact and in analytical closed forms. They are usually 
used as basic solutions of conformal transformation methods as we will dis- 
cuss in the next section. This method will be used in Chapter 5 to calculate 
the flow field behind a circular flameholder. The mathematical formulation 
of this method will be outlined in detail there. 

2.1. 1.2. Conformal Transformation 

For the flow in a duct, around a corner or past an elliptic cylinder, the 
inviscid flow solution can no longer be simply obtained by the method of 
images. By expressing two-dimensional coordinate in the form of complex 
variables, conformal transformation can then be employed to transform a 
given geometry into an upper half plane or a circle for which solutions are 
already known. Several well-known examples in this category are Schwarz- 
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Christoffel transformation which transforms a closed polygon into an upper 
half plane, Joukowski transformation (Milne-Thomson, 1971) which 
transforms an airfoil or an elliptic cylinder into a circle and Theodorsen's 
method (Theodorson, 1932) which transforms a single arbitrary-shaped body 
into a circle. Among these methods, Schwarz-Christoffel transformation is 
especially useful in our numerical model because it can be easily used to 
apply to a variety of combustion configurations such as a sudden-expansion 
duct, a backward-facing step and a cavity. The method transforms a given 
polygon in the physical plane (z-plane) into an upper half plane (f -plane) 
by the following differential equation. 

Fit) = = K{ t-a)" _l ( t-b) • ( f-c) - _1 • • • (2.10) 

where a , 0 . y ... are interior angles of a simple closed polygon of n vertices 
so that <* + /? + 7 +... = (n-2) n. 

the parameters a,b,c a < b < c, are the coordinates of n points on 

the real axis in the t -plane corresponding to n vertices in the z- 
plane. 

K is a constant to be determined ; it may be a complex number. 

The relation between z and t can be found by integrating equation 
(2.10). The result usually takes an integral form and has to be evaluated by 
numerical integration except for some particular geometries (for example, a 
parallel duct and a splitter plate) for which analytical formula can be 
obtained. The velocity field in the physical plane is evaluated by first calcu- 
lating in the transformed plane and then multiplying the transformation 
function obtained above. 
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This method is easy to apply. In some cases, however, integration con- 
stants become very difficult to find and the numerical integration is difficult 
to perform accurately (Sand, 1984). For some specific geometries such as a 
double-sided sudden expansion duct used in a dump combustor, the 
transformation function has the form of a strong exponential increase or 
decay regarding channel length (Giovaninni, 1984 ). This restricts the 
dimension of the flow domain which we are interested in due to the limita- 
tion of computer accuracy. 

When conformal transformation is used to find the inviscid flow solution 
with vortex method, it can reserve the grid-free principle of the random vor- 
tex method and is thus devoid of numerical diflusivity. But the deformation 
of the shape of vortex blobs due to conformal mapping is speculated to 
create certain amount of inaccuracies. 

2. 1.1. 3. Integral-Equation Methods 

The previous two methods can only be applied to two-dimensional flow 
about a single body. For flow about a three-dimensional body or two- 
dimensional multiple bodies, integral-equation method seems to be an excel- 
lent way to And the potential flow solution. The ideal of this method is to 
reduce the problem to an integral equation by a singularity (source-density 
or vortex- density) distribution on the surface of the body. According to 
Hess Sc Smith (1987), the resulting integral equation for the source-density 
distribution a (p) is 

2ttj(p) )"(*)«*» = — n (?) ' U m (2.11) 

where n denotes a derivative in the direction normal to the boundary 
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surfaces at point p, q is a point outside the body and is the velocity at 
infinity. The conditions under which a solution of the above equation can be 
obtained are general. The shape of the body does not need to be analytic 
and may consist of several disjoint surfaces. Detailed solution procedures 
and theories of two- or three-dimensional flows can be found in the works by 
Kellogg (1929) or Hess & Smith (1967). 

Since this method is used to solve for the potential flow only, the total 
inviscid velocity u in equation (2.7) has to be decomposed into a rotational 
velocity component u y and a potential velocity component Up by the well- 
known orthogonal decomposition, 

u = Uy + Uj, (2.12) 

where 

Up = V 

and 

V x = V z * v = - ( (2.13) 

V X Up = V 8 'frp = 0 (2.14) 

the boundary condition becomes 

u n = (iip + Uy) n - 0 (2.15) 

or 

Up ■ n = -Uy n (2.16) 

The rotational velocity component u y can simply be found by summing up all 
velocities induced by vortex blobs according to equation (2.9) while the 
potential velocity component Up can be evaluated by applying equation 
(2.11). The method was used by Chorin (1973) to find the potential velocity 
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in his study of flows past a circular cylinder and it will also be used in our 
future studies. 

2.1. 1.4. Finite -Difference Methods 

When the flow domain has a more complicated geometrical form, finite 
difference method offers a most straightforward way to approximate partial 
differentiation operators to solve the problem. Equation (2.8) can be 
directly discretized by using central-difference to approximate the Laplace 
operator 


£F< *i-U + *<+i.j + *<j - 1 + *U-i “ 4 *ij ) = £ ij (2.17) 

where i. j are indices of meshes and h is the mesh size. 

2.1. 1.4.1. The Cloud-in-Cell Method 

In order to retain the Lagrangian description of the vorticity field but 
solve the Poisson equation for the velocity field on a Eulerian mesh, a 
method called cloud-in-cell was used by Christiansen (1973) and Baker 
(1979) in research on the interaction of regions of vorticity and the roll-up 
of vortex-sheet. This method uses area weighting scheme to redistribute 
vorticities inside a given mesh cell to the four mesh points at the corners 
and solve equation (2.17) by fast Poisson solvers. The velocity of each vor- 
tex blob is then determined by central-difference differentiation of stream 
function and bilinear interpolation to convect itself in the Lagrangian coor- 
dinate. The major advantage of this method is the reduction of computa- 
tional time from 0(N y ) of direct calculation to G(N y ) + O(MlogM), where N y 
is the number of vortex elements and M is the number of grid points. This is 
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especially useful when the number of vortex elements reaches a very high 
value. However, the method introduces some additional numerical errors 
due to the anisotropic nature of the vorticity distribution and interpolation. 
The angular moment of vorticity also does not conserve. For detailed dis- 
cussion of this method and error considerations, see papers by Christiansen 
(1973) and Baker (1979) as well as the review by Leonard (1980). 


2.1. 1.4.2. The Successive Ovei^ Relaxation Method 


There are some cases for which good accuracy is desired and the 
number of vortex elements is not too large. To prevent the anisotropic 
nature of the vorticity distribution, one can solve only potential flow com- 
ponent by the finite-difference method and calculates the rotational velocity 
component by direct vortex interactions. The equation to be solved for 
potential velocity is 


V 2 i p = 0 

with boundary condition 


3Sp , 

-q£~ I W - Uj, n | u, 


“T*vn 


(2.18) 


(2.19) 


or they can be expressed in terms of stream function as 

V 2 *, = 0 (2.20) 

and 


Qg lw = = ' u vn (2.21) 

where n is the unit vector normal to the wall while s is the unit vector along 
streamline and u^ is the velocity induced by vortex blobs in the normal 
direction. 
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By using five-point approximation for Laplace operator, the finite- 
difference form of equation (2.18) and (2.20) can be obtained in the same 
form as equation (2.17) except that the right hand side of the equation 
equals to zero instead of . Equation (2.18) and (2.20) can be solved by 
either direct methods ( e.g. Fourier method (Dorr, 1970), EVP method 
(Roache, 1971)) or iterative methods ( e.g. SOR method (Young, 1954) and 
ADI method (Peaceman & Rachford, 1955)). Among these methods, the SOR 
(successive over-relaxation) method seems to be the easiest one to under- 
stand and to program with compatible efficiency for complicated flow 
domains. It is thus chosen here to solve the Laplace equation for the poten- 
tial velocity. 

By viewing equation (2.18) and (2.20), the potential velocity Up can be 
solved either from velocity potential or from stream function. Since equa- 
tion (2.21) will give us the Dirichlet boundary condition which is easier to 
deal with numerically, we would determine to solve the Laplace equation in 
terms of stream function. According to the successive over-relaxation 
method, equation (2.20) can then be written in the discretized form as 

+ *rp i + + *,Vi - 4 *u] (2.22) 

where superscript n, n+1, refer to number of iterations and u is the over- 
relaxation parameter 

The boundary condition (equation (2.21)) can also be written in the 
finite-difference form as 


As 


-'Uvr, 


(2.23) 


or 
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= ** -Uvn ■ AS 


(2.24) 


where subscript k and k+1 are indices of successive grid points along solid 
walls and A s is the spatial increment between two points. 

The iteration of the SOR method will converge only when the value of 
the over-relaxation parameter u is between 0.0 and 2.0. In a complicated 
region, there is no analytic formula for predicting the optimal value of o. It 
thus has to be determined by numerical experiments. However, the analytic 
formula for the problem in a rectangular domain of I x J grid points derived 
by Frankel (1950), 


u 0 


2 


1 - VT^S 
X 


(2.25) 


where 


X 


cos( ) + COS ( J ~ i ) 


(2.28) 


is used to estimate the optional parameter u for the initial guess. 

The formula of equation (2.22) can be used to calculate stream func- 
tions for the interior grid points in the flow domain as long as the shape of 
the mesh remains a square. In the case of a sloping wall or a curve boun- 
dary, grid cells near the boundary become irregular-shaped (see Figure 1), 
the finite-difference approximation of Laplace operator has to be modified in 
order to handle the irregularity of grid cells. A second-order Taylor series 
expansion is utilized in order to approximate the partial differentiation 
operator as : 
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(2.27) 


dH 

2 

*t-U 

*R *U 

dx 2 

h 2 

1 + 4r 

#r( 1 + &R ) &R 


dH _ 2_ *B 

dy 2 h 2 1 + $b &b(1 + "A B ) 



(2.28) 


thus 


V 2 4- 


dH_ 

dx 2 + 9y a 


2 

/i a 




*R 


1 + 1>J» ‘tfj? ( 1 ) 


^iiL. + ( J_ + -L. ) * 

+ (1 + &B ) ^J? 


(2.29) 


where = I CR I /h 
= I CB I /h 

'I'g Is the stream function at point B calculated by equation (2.24) 
4^ is the stream function at point R calculated by interpolation. 


The discretized Laplace equation used in the SOR method for irregular 
grids becomes then : 


♦.v 1 = 


+ 


u 






*b 


♦i 


n 


& r ( 1 + “So ) 1 


* l 


1 + ^R &b( 1 + g ) 


(2.30) 


In order to calculate the stream function at each grid point on the 
boundary, the velocity induced by vortex blobs normal to the wall u^ needs 
to be evaluated. When the boundary is parallel to x- or y- axis, the normal 
velocity component is simply the y or x component of the total induced velo- 
city. In the case of a sloping wall or a curve boundary, it can then be simply 
found by coordinate transformation as 
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u sintf + v costf 


(2.31) 


«n = 

where is the angle measured counter-clockwise from the x-axis to the 
tangent of the boundary at certain grid point, and u, v are velocity com- 
ponents in the x and y direction respectively. 

For internal flows, besides boundary conditions at walls, inlet and 
outflow condition, also have to be specified in order to be able to solve the 
Laplace equation. The inlet flow condition is usually given and its value can 
be directly assigned to inlet grid points. But outlet flow condition is often 
not known o prior. It has to be evaluated with the entire flow field at the 
same time. Thus, instead of using a prescribed velocity profile, a less res- 
trictive condition d'I'/flx = Ois suggested here as the boundary condition 
at the exit of the duct. This condition has been reported (Paris & Whitaker, 
1965), and is also found in our study to be quite successful in obtaining 
accurate results, as we will discuss in Chapter 4. 

For external flow such as would occur around a blunt body, the flow 
domain often extends to infinity. In this case boundary conditions at the 
upper and lower boundaries of the computational domain have to be 
specified. Since the mesh size used in the potential flow calculation does 
not have to be very small, one can always specify an undisturbed velocity 
condition at upper and lower boundaries far enough from the body. For a 
parallel flow past a blunt body we can use U = U , where U is the velocity 
upstream , or constant stream functions as boundary conditions. 

Once all boundary conditions have been specified, one can start to 
iterate either equation (2.22) for square meshes or equation (2.30) for non- 
square meshes by successive over- relaxation method until Max I ' 

’fry” I < t .where r is the desire accuracy and I, J are the number of grid 
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points in x- and y-direction respectively. The number of iterations required 
for each time step will depend on the specific geometry under consideration. 
For the backward-facing step with a contraction inlet, as will be studied in 
Chapter 4, it ranges from 20 to 40 depending on velocity fluctuations of each 
time step. 

The potential velocity component of each vortex blob may be found by 
differentiating stream function with central difference and by interpolating 
velocities at grid points at four corners of a mesh using area weighting 
method (Figure 2) 


= 



(2.32) 


The total convection velocity of each vortex blob is then obtained by adding 
rotational velocity u y to the potential velocity Up according to equation 
( 2 . 12 ). 


2.1.2. Viscous Effects 

In addition to the zero normal velocity boundary condition of an Euler 
flow, the tangential velocity component also has to vanish for viscous flow 
due to finite value of fluid viscosity. The no-slip boundary condition is 
satisfied by creating vorticities at rigid walls to cancel non-zero tangential 
velocities. The diffusion equation 


§f = Re' 1 V 2 t (2.33) 

is present here together with the Euler equation in order to constitute the 
complete vorticity transport equation. These two major components of 
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viscous effects in our model are now briefly discussed as follows : 

2. 1.2.1. Diffusion Equation and Random Walk 

The diffusion process of vorticity in viscous flow is characterized by the 
diffusion equation with Re'* as the coefficient. The solution of diffusion 
equation in free space is a Gaussian distribution of mean zero and variance 
2t / R 


G(x,t) = ^/^exp(-^-i) (2.34) 

For conventional time-dependent flow calculation, the time interval (0,t) is 
always split into n small time steps, each with interval At = t/ n. Since the 
vorticity field in our model is represented by a number of vortex blobs, we 
can let them undergo random walks which are independent Gaussian ran- 
dom variables with mean zero and variance 2 A t / Re. This is because two- 
dimensional flow random displacements in two directions are also indepen- 
dent to each other. The final displacement after n time steps thus also has 
a Gaussian distribution with mean zero and variance 2t / Re. Thus the dis- 
tribution of vortex blobs at time t is given by equation (2.34) and it satisfies 
the diffusion equation. The location of each vortex blob at next time step 
can be found by 


^ dt + 

(2.35) 

V r* 1 = Vf + / Vi dt + rjty 

(2.38) 


where rj- lx , rj . are Gaussian distributed variables with zero mean and vari- 
ance 2 A t / R. 
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Since the computational time required for velocity calculation at cer- 

p 

tain location is proportional to N , where N is the number of vortex blobs, a 
higher order time integration usually requires the evaluation of velocities at 
certain intermediate time steps and thus increases computational time 
dramatically. It is certainly not worth pursueing if higher accuracy is not 
critical. First order Euler's time integration is often used to carry out the 
displacement calculation of vortex blobs. 

2.1. 2.2. Vorticity Generation at Boundaries 

In Chorin’s paper (1973) of random vortex method, the tangential velo- 
city at boundary is eliminated by creating vortex blobs of total vorticity ( = 
u w h at each equally spaced check point along the wall, where u w is the 
tangential velocity at wall and h is the space between check points. By 
doing this, the rate of convergence near boundary is slow and the computa- 
tion effort is substantial especially in interior flow problems. Therefore, a 
vortex sheet approximation of boundary layer was developed (Chorin, 1978) 
to overcome these problems. 

According to boundary layer approximation, the velocity gradient in 
streamwise direction is small compared with that in the normal direction. 
The vorticity equation is thus simplified as 

If + ( U -V) { = Re-0. (2.37) 

where 

( = -% (2.38) 

The diffusion component of vorticity equation becomes one-dimensional, 
therefore, only random walk m y-direction is needed to simulate the 
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diffusion process of vorticity. Detailed formulations of convection velocity 
inside boundary layer resulting from above equations have been derived by 
Chorin (1978). Vortex sheets are designed to induce velocity field only below 
themselves. Each time step new vortex sheets are generated at walls 
according to amount of vorticities needed to satisfy the no-slip condition. 
They are then subject to random displacements of zero mean and variance 2 
A t / Re in y-direction. In order to attain certain computational accuracy, a 
maximum vorticity T „ allowed is selected for vortex sheets. This results 
in generating several vortex sheets at certain check points along walls. 

Since boundary layer approximations are valid only inside the boundary 
layer of which the thickness is of the order of 0(Re~^), vortex sheets approx- 
imation can be applied in that small region only. The numerical boundary- 
layer thickness <5 is chosen to be a multiple of standard deviation (2 A t / 
Re)^ in order to minimize the loss of vortex blobs across boundaries (Chorin, 
1978). When vortex sheets move outside the boundary-layer, they are 
turned into vortex blobs and induce rotational flow outside. On the other 
hand, when vortex blobs move into boundary-layer, they are changed to vor- 
tex sheets and follow the motion of boundary-layer flow. As we stated above, 
vortex sheets only induce velocities below themselves. Therefore, when cal- 
culating velocities at certain points, we don't have to include the effects 
induced by vortex sheets unless the point of interest is inside the numerical 
boundary layer. The computational effort is thus tremendously reduced 
especially in internal flows for which the number of vortex sheets is large. 
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22. Turbulent Reacting Flows 

As a consequence of all the complexities associated with turbulent flow 
in the presence of chemical reactions, the analysis must be based on an 
array of simplifying assumptions. In particular, we are interested here in 
the case of a perfectly premixed reacting medium where chemistry is so fast 
that combustion is controlled primarily by turbulent transport. As already 
stated earlier, flame front is modeled as an interface between premixed 
reactants and burned products that is significantly thiner than smallest tur- 
bulent length scale. Thus the only effect of turbulence on the flame is mani- 
fested in wrinkling its front while its integral features remain unchanged. 
According to Ghoniem et al(l982), under such circumstances, the flame 
motion is described as follows 


and 


V • u = t(jf ) 


(2.39) 


-fif- = n + u (2.40) 

where r^ is the position vector of the flame front, S u is normal burning speed 
assumed to be a prescribed constant and n is the unit vector normal to 
flame front, while 

t ( Tf ) m u, -5 ( r -T f ) 

u. = j(v- 1) Su 

models the local rate of expansion by volume sources. In the about V = v^ / 
v u is the specific volume ratio across the flame front. 
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2.2.1. Flame Ftont Location : the SLIC Method 


The non-reacting flow is described in the Lagrangian coordinate by fol- 
lowing the motions of vortex blobs, while the flame front is traced at the 
same time in every time step. 

The shape of flame front in a turbulent flow always exhibits such a com- 
plicated characteristics that it is difficult to represent it by a single interpo- 
lation function, especially when cusps occur. Quite often, one finds also the 
appearance of the so-called unburned ’’islands" and "pockets" (Bray, et al, 
1984) behind the flame front due to the entrainment of the unreacted 
medium. Therefore the flame, which is not even simple-connected becomes 
more difficult to describe by an analytical function. 

The method used here to determine the location of the flame front is 
based on the SLIC(Simple Line Interface Calculation) algorithm (Noh & Wil- 
liam, 1976 and Chorin ,1978 ). According to this the flow domain is first 
divided into a number of square ceils. Volume fraction of burnt medium 
inside each cell f is calculated and the geometry of the interface is deduced 
from the value of f in the cell and its four neighbors. Interface geometries 
inside cells of partially burned are classified into four cases : (i) vertical 
interface, (ii) horizontal interface ,(iii) corner and (iv) neck, as illustrated in 
Figure 3. 

2.2.2. Self -Advancement and Huygens Principle 

Since no explicit function is given to describe the shape of flame front, 
the normal direction of every point along the flame front cannot be simply 
determined. The self-advancement of the flame in its normal direction is 
computed by an implementation of the Huygens principle. Interface inside 
the cell is moved at a constant flame speed into eight different directions, 
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a 

each 45 apart. The maximum value among eight new f numbers is chosen 
to be the resulting volume fraction. New values of volume fraction f are 
used to determine the interface geometries and locations. The change of 
volume fraction f inside each cell is also recorded to evaluate the strength 
and core radius of source blobs. 

An example of flame propagation in its normal direction can be easily 
obtained by igniting premixed reactants inside a rectangular duct at the 
center of its end wall. The shape of flame front remains a semi-circle until it 
reaches side walls. It then becomes flatter and flatter due to the increase of 
radius of curvature. In Figure 4, we display the numerical simulation of this 
case for every 20 time steps. Numerical parameters used are h = 0.05, A t = 
0.1 and S u = 0.25, where h is the cell size, A t is time step, and S u is the 
prescribed laminar flame speed. 


2.2.3. Velocity Field of Reacting Flows 

The dynamic effect of exothermicity during the combustion process was 
properly modeled by the use of volumetric source blobs located along flame 
front. The strength and core radius of each source blob is evaluated 
according to the change of volume fraction f in the cell by normal burning. 
The velocity field induced by source blobs in free space is described by 


u . - £ A, 

Z — Zj 

1 

• 2n / 

i=i 77 max( 

Z - Zj 

,r. ) (* " ) 


where Aj is the strength of the source blob 

r g is the core radius of the source blob. 


(2.41) 
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A source blob is in concept analogy to a vortex blob, or in terms of com- 
plex function, a source blob is a vortex blob with imaginary strength. Equa- 
tion (2.41) has the same form as equation (2.9), except that the velocity field 
produced by source blobs is potential and curl free, while that induced by 
vortex blobs is rotational, but divergence free. 

If the method of image or conformal transformation can be used to find 
the solution of the flow field, a similar way used in section 2.1. 1.1 and 2.1. 1.2 
for vortex blobs can also be applied to source blobs with the exception that 
image source blobs do not change signs of their strengths. When integral- 
equation or finite-difference method is applicable, the velocities induced by 
source blobs have to be added to the total velocity in equation (2.12) as 


U = It* + Up + u. 

(2.42) 

Vxu*=V 8 * v = -£ & V ■ u* = 0 

(2.43) 

V x Up = V 8 = 0 * 7 ■ Up = 0 

(2.44) 

Vxu, = V 8 = 0 & 7 u, = e 

(2.45) 


The boundary condition for solving potential velocity becomes 


or 


it • n = (up + it* + it, ) • n = 0 

Up n = - (u„ + u,) - n 
= — (^^vn + ' u »n) 


(2.48) 


(2.47) 


The potential flow component can then be solved by methods outlined in 
section 2. 1.1.3 and 2.1. 1.4 . Total velocities of reacting flows can be obtained 
by adding up u g in equation (2.41) and u y in equation (2.9) to Up. The 
overall procedures for calculating reacting flows are depicted in Figure 5 by 
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block diagrams. 

Flame propagation is usually initialized by a point ignition or a line igni- 
tion. Sometimes a region with completely burned products is introduced in 
order to save the computational time of establishing a continuous front. 
Initial interfaces inside cells are moved by the normal burning speed and 
local convection velocity according to equation (2.40) in two fractional 
steps. The number of cells used in the reacting flow calculation may be 
large. In order to reduce the computational effort, flow velocities are calcu- 
lated only when either cells themselves or at least one of their eight neigh- 
boring cells are active. A cell is considered to be active at the specific time 
step when the change of volume fraction of burned products in it during the 
normal burning process is finite. This algorithm lets us move flame front 
every time step correctly, without going through all cells in the computa- 
tional domain. 
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Chapter 3 


Rudiments of Flame Aerodynamics in Premixed Gases 


Non-steady flame propagation in a variety of combustion configurations 
have exhibited many interesting features which include the coherent struc- 
tures in reacting flows, the breakdown of laminar flames and cellular struc- 
tures of flame front stabilized by a bluff body. Each of these phenomena 
involves a large set of complicated mechanisms which are difficult to deal 
with, even individually. 

In premixed combustion fluid mechanics plays an important role in 
determining the reaction rate, flame shape and other combustion processes. 
By relating the salient features of flame propagation to the aerodynamic 
properties of the flow field, we are able to deduce three fundamental 
mechanisms which physically are capable of interpreting most flame 
phenomena in premixed gases, and mathematically are particularly con- 
nected with our numerical model. 

While the detailed calculations of turbulent flow fields and the develop- 
ment of flame front in both confined and unconfined regions will be carried 
out in Chapter 4 and 5, we will, first, in this chapter use an ideal potential 
flow model to illustrate the individual or combined effects of these funda- 
mental mechanisms on the flame propagation in simple geometries. Exten- 
sions are possible to cover most aerodynamic features in premixed gases. 
The formation of tulip shape flame in a closed duct serves as another exam- 
ple of our illustrations. 
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3.1. Three Fundamental Mechanisms 


The flow fields of most practical combustion systems, whether turbulent 
or laminar, often consist of rotational flows of various length scales. The 
rotational flow field which may result from wall effects or density gradient 
across flame front can always be represented in the form of vortices. These 
small or large vortex structures in reacting flows have the effect of entrain- 
ing more fresh reactants into to the reaction zone and then modifying the 
mixing and reaction rate of combustion process. The rotary motion of 
eddies thus plays an important role in the determination of flame shapes. 

Combustion is a process associated with the transformation from reac- 
tants into products at a certain reaction rate. The normal burning velocity 
of flame front has a strong dependency on local properties of reacting gases 
as well as the structure of the flame front. Accordingly, no unique function 
has yet been found to describe the relation between the normal burning 
velocity and its dependent variables with satisfactory results. However, 
models which assume constant burning velocity are found to be very suc- 
cessful in describing many flame propagation phenomena. The effect due to 
constant normal burning is also considered here to be one of the basic 
mechanisms in the investigation. 

The dynamic effect of thermal expansion result from chemical reaction 
is modeled by volumetric sources located along flame front as described in 
previous sections. These volumetric sources not only affect the flow field 
but also modify the shape of flame front due to the interaction between 
them.. 

The three fundamental mechanisms of flame propagation -- (i) rotary 
motion of vortex structures; (ii) self-advancement of flame front at a 
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prescribed constant normal burning speed; and (iii) the concomitant action 
of volumetric sources along flame front used to model the thermal expan- 
sion, have already been brought up by Oppenheim & Ghoniem (1983). We 
would like to use these mechanisms to investigate some aerodynamic 
features of flame propagation in premixed gases in more detail. 

In this section the simple case of a straight flame front normal to the 
walls of a parallel duct is used to study the deformation of the flame front 
due to the individual or combined effects of three fundamental mechanisms 
stated above. To make the problem even simpler, we assume the flow field is 
potential, i.e. 


where 


V z $ = 0 


(3.1) 


with boundary condition 


u = V • $ 


(3.8) 


8 £ 

dn 


0 


(3.3) 


The velocity field induced by a potential vortex located at the center of 
the duct is formulated by the Schwarz- Christoflel transformation. The 
transformation function of a parallel duct is found to be 


r«) = ^ = n( (3.4) 

The relation between z and f can be easily obtained by integrating the above 
equation as 
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or 


= exp ( rr * ) 


(3.5) 


z = ~logf (3.0) 

For a potential vortex with strength T located at z Q in physical plane, the 
complex velocity in the transformed plane is specified by 


where 


*(<) 


iT_ 1 1 

2tr f " <o < + 

<rg. 

*((* + <$ ) 


(3.7) 


= exp ( nz 0 ) (3.8) 

The velocity field in the physical plane can then be obtained by taking the 
conjugate of the product of the transformation function and the complex 
velocity as 


u(z) = W\z) = 


(3.9) 


where ( )• denotes the conjugate of a complex number. To perform the com- 
putations in terms of nondimensional parameters, we choose the coordinate 
system shown in Figure 8. a. The height of the duct H is chosen to be 1.0, 
thus the coordinate of the vortex is z Q = 1/2 , The strength of the vortex T is 
set to be 2.0 so that velocity at the nearest wall points induced by the 
potential vortex is equal to unity. The normal burning velocity is taken as 


3, = 0.4 

while the density across the flame front is 
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v fc 

V = — = 3.0 

▼u 

where subscript u and b denotes the unburned and burned gas respectively. 
The source strength per unit length of the flame front is 

u, = j(V -1)5* = 0.4 (3.10) 

and the time step is taken to be 

Af = 1.7857 x 10' 3 

in order to satisfy Courant condition in the computation of flame propaga- 
tion. 

Figure 7 displays the deformed flame fronts and velocity vectors on 
both unburned and burned sides after 80 computational time steps. The 
scale of velocity vectors is shown at the right corner below the first plot. 
The dash line represents the initial flame front while solid lines show the 
contours of deformed flame fronts. 

Figure 7.a applies to the case of the distorted flame front due to the 
rotary motion of the vortex alone. A large deformation of the flame shape as 
well as an increase of flame area resulting from the interaction of the flame 
with the vortex element can be observed from the figure. When the normal 
burning speed is added to the motion of the flame front, the velocity fields 
on both sides of flame front of course remain unchanged, while the flame 
front advances in its normal direction, smoothes out the wrinkles and forms 
a cusp as the result shown in Figure 7.b. 

If the flame front in Figure 7.b is compared with that of the initial 
straight interface propagates at the same flame speed without the vortex 
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motion, one finds that the burning rate, which is defined here as the number 
of burned ceils consumed per time step, of the former case is about two 
times higher than the latter. The increase in the flame area due to the 
stretch of vortex motion obviously augments the burning rate. That tur- 
bulent burning velocity is larger than the laminar one due to the wrinkling 
of flame front by the interactions of turbulent eddies has been observed in 
most turbulent combustion. 

Figure 7.c displays the flame front due to the combining effects of all 
three fundamental mechanisms. The basic flame shape remains unchanged 
from Figure 7.b except that flame front advances more into the unburned 
gas by the augmented flow field. Velocity vectors on the unburned side are 
found to have significant changes due to effects of volume sources. The 
directions of the velocity vectors also indicate the divergence of streamlines 
outside flame fronts. 

While the velocity vectors on the burned side are found to be 
unchanged by the action of distributed volume sources. This is because, as 
interpreted by Oppenheim & Ghoniem (1983), that physically the domain of 
burned gases contains a point of zero translation velocity, and mathemati- 
cally the domain of burned gases in the transformed plane is closed by the 
flame front, continuous volume sources distributed along the front thus do 
not change the flow field inside. 

To illustrate this interpretation, the computation of flame front and 
velocity vectors with same conditions in the transformed plane are per- 
formed. Figure 8 is the corresponding results of Figure 7 in the f - plane. 
The enlargement of the inner portions of Figure 8 is shown in Figure 9. The 
scale of velocity vectors of each figure is also shown at the right hand 
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corner below the first plot. The closed feature of the burned gases can be 
clearly noticed in the figure. 

Although the calculation is done in a confined region, similar results as 
those shown in Figure 6 - 8 in an unconfined domain can also be obtained 
except that the zero normal velocity boundary condition is excluded. Thus 
the flame front does not necessarily remain a right angle with duct walls. 
Cusps will still exist at the flame front and form singularities in the flow 
field. The result is similar to the cellular structures of confined or 
unconfined flames stabilized by a bluff-body holder. A phenomenological 
sketch by Oppenheim (1983) (Figure 10) is also found to be simulated well by 
our numerical model. 


3.2. Flame Propagation in a Duct with Closed End 

The phenomena of flame propagation in a tube have been investigated 
for over a century (Mallard & Le Chatelier, 1883, Maxworthy, 1962, Zel’dovich 
1980). Most of the experimental or theoretical works have concentrated on 
the steady plane flame and result in very complete mathematical theories 
(Sivashinsky, 1983) in this area. The phenomena of non-steady flame propa- 
gation have also been studied by many researchers; however, results are 
less satisfactory. 

One of the interesting features of non-steady flame propagation in a 
closed duct is the formation of the so-called "tulip shape" flame. If the 
length of the duct is not too long, it has been found (Ellis. 1928, Guenocho 
1964, Steinet et al, 1982 and Dun-Rankin et al, 1984) that, when the 
premixed reactants inside a closed duct are ignited at one end, the flame 
initially exhibits a semi-cylindrical shape and then flattens into a planar 
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front after it reaches side walls. As the planar front propagates further 
downstream with a noticeable decrease in burning rate, cusps begin to form 
along the front. The center cusp of the front then becomes more and more 
pronounced and finally develops to a "tulip” shape flame as the flame 
approaches the closed end. 

Since no quantitative results have been obtained on this problem, the 
causes of tulip shape flame remain controversial at this time. Flame stabil- 
ity, pressure wave and flame driven flow field are three most prevailing 
speculations about this interesting phenomenon. 

The measured maximum flow velocity induced by the flame in a closed 
channel is found to be of the order of several meters per second. The Mach 
number in this case is small, thus pressure wave will not have a significant 
effect on either the flow field or the flame front. For stoichiometric 
methane/air mixture, the Lewis number, which isflned as Le = D /k, approxi- 
mately equals 1.0. The value is near the thermo-diffusive instability thres- 
hold of flame fronts, hydrodynamic instability therefore becomes a dom- 
inant factor. According to Markstein's theory, the flame is unstable to 
long-wave disturbances but is stable to short-wave ones. The dimensions of 
ducts or tubes used in experiments are always small. The planar flame front 
bounded by narrow walls can thus prevent the appearance of long-wave dis- 
turbances and becomes hydrodynamically stable. 

We therefore believe that the flow field driven by the flame is the major 
cause of the "tulip" shape flame. As experimental photography showed, at 
the initially stage of flame propagation the effect of flame driven flow is 
small. The flame behaves like a free-propagating front until it "feels" the 
existence of the end wall and then flattens. This early part of flame propa- 
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gation can be easily explained or modeled by the Huygens Principal as we 
discussed in last section. 

The transitional flame front is modeled here by the use of a simple con- 
ceptual model. It considers a planar flame front propagating in a duct 
toward the closed end wall. The coordinate system used for the computa- 
tion is shown in Figure 6.b. In view of the symmetry of the flame, the compu- 
tation is performed only on the upper half of the duct. According to 
ZeTdovich's analysis (i960), the flow field ahead of the flame is potential 
while that on the burned side is rotational. The inviscid flow model used in 
the last section is also chosen here to calculate the flow field. The rota- 
tional flow part is simulated by creating a large scale vortex behind the 
flame. 

To calculate the velocity field in the duct, the Schwarz-Christoflel 
transformation is used to map the computational domain into an upper half 
plane. The transformation function for the geometry under consideration is 
found to be 


HO = ~ (3.11) 

After integration, we have the relation between z and f , as 

(■ = cosh[rr(i -z)] (3.12) 

Assumed the location of the vortex at the physical plane is z Q , the complex 
velocity in the transformed plane can be found by the formula in equation 
(3.7) with = cosh[rr (i-z )]. The velocity field in the physical plane is then 
obtained by 
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(3.13) 


u(z) 


w\z) = (*({) /’({)]' 
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For the sake of convenience, tbe width of the duct H is also chosen to be the 
reference length. In order to allow the planar front to have enough distance 
to develop into a tulip shape, the flame is initially located at 2.5 H away from 
the end wall. The vortex is then located at the center of the half duct and 
0.5 H behind the flame front. The strength of the vortex is chosen to be -2.0 
so that the velocity at the nearest wall point is equal to unity. For con- 
sistency, the values of laminar flame speed and density ratio are chosen to 
be the same as those used in the last section, i.e. S u = 0.4 and V = 3.0 

If the location of the vortex is Axed in space with time, the effect of 

rotary motion acts on the flame front will be gradually reduced as the flame 

keeps on propagating downstream. In order to keep up the effect of the vor- 

-3 

tex, a constant velocity with the magnitude 7.5 x 10 to the right is chosen 
to impose on the motion of the vortex after a number of numerical experi- 
ments. Computational time step A t is taken to be 1.7857 x 10 to satisfy 
the stability condition. 

Figure 11 displays the development of flame fronts and locations of vor- 
tices at the 40th, 120th and 200th time step. The initial flame front is 
represented by a dash line while solid lines delineate tbe contour of flame 
fronts at subsequent time steps. Figure 11. a, b k c. are results of three 
cases corresponding to Figure 7. a, b & c, i.e. the effect(s) due to (a) the 
rotary motion of vortex alone; (b) rotary motion and self- advancement in 
normal direction, and (c) all three fundamental mechanisms. 

By the action of the rotary motion alone , flame fronts are stretched 
and exhibit peculiar shapes. The added flame advancement in the normal 
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direction smoothes out the front and forms flame shapes similar to the tulip 
flame observed in experiments. When the effect of thermal expansion is 
modeled, the axial motion of flame front is found somehow retarded and the 
fold of the front is deeper into the burned gases. This is because that the 
domain of unburned gases surrounded by the flame front and duct walls is 
now closed. The flow field induced by distributed sources along flame front 
thus only affect the velocity field on the burned gases without modification 
of that on the unburned gases. 

Velocity vectors and flame fronts at the 200th time step are shown in 
Figure 12, while Figure 13 shows the corresponding plot in the transformed 
plane. 

Although the model used in this section is only phenomenological and 
over-simplified, we believe that we have already demonstrated some impor- 
tant features of the unsteady flame propagation in a closed duct by such a 
simple model. A more detailed numerical model which calculates the flow 
field induced directly by the changes of pressure and density from the 
exothermicity of chemical reaction in a closed domain is currently under- 
way. The location of flame front is determined by following the flow velocity 
produced by itself. Meanwhile, the experimental measurements of the 
detailed flow field by the LDA technique is also in process. We hope all of 
these efforts will bring us to a better understanding of the formation 
mechanism of tulip shape flame. 
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Chapter 4 


Confined Flame: 

A Turbulent Flame Stabilized Behind a Rearward-Facing Step 

Combustion of premixed gases in the wake of a rearward-facing step 
has recently attracted deservedly broad attention, as manifested by experi- 
mental investigations of Ganji & Sawyer (1980), Pitz Jc Daily (1983), 
Shepherd, et al (1982) and El-Benhawy, et al (1983), as well as by the numer- 
ical studies of Ghoniem, et al (1982) and Ashurst (1981). The reason for this 
lies in the following unique interesting features of the system: 

(1) it provides proper means for mass and heat recirculation, a process of 
particular advantage to lean combustion, enhancing thermal efficiency 
and reducing pollutant emission; 

(2) it promotes strong interaction between the turbulent flow field and the 
combustion process, forming a closed loop feedback mechanism; 

(3) it exploits the intrinsic flow instability to stabilize the flame; 

(4) it enhances mixing, leading to the establishment of close coupling 
between the combustion process and flow instabilities. 

However .because of the same reason, the instability of the flame is of 
crucial important 3 ince it is so dependent on various operating conditions, 
such as equivalence ratio, inlet flow velocity, pressure, and temperature. As 
the flame is stabilized by the recirculation, bringing hot products into close 
contact with the reaction zone furnishing thereby a continuous supply of 
heat for ignition of fresh mixture, the turbulent field plays a major role in 
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determining the geometry of the flame front. On the other hand, the expan- 
sion associated with the conversion of reactants into products acts to 
modify the turbulent field, e.g. the recirculation (reattachment) zone length 
is reduced to almost one half of its value in the non-reacting flow due to gas 
acceleration across the flame front. Meanwhile, the free-stream flow and 
the recirculation vortex interact with the set of trailing eddies of the mixing 
layer dominated by large scale structure, enhancing the interaction 
between combustion products and reactants and acting locally as a heat 
exchanger that augments the temperature of the reactants. There is a wide 
range of frequencies detectable in the energy spectrum of the mixing layer 
corresponding to the size and motion of these eddies. However, some fre- 
quencies can become dominant if the mixture quality is adjusted, resulting 
in flow instabilities that cause the flame to oscillate (Keller, et al, 1982). 
This process emphasizes the close coupling between fluid mechanics and 
combustion in the system under study, as well as the major role of the vorti- 
cal flow structure that has to be, therefore, properly taken into account by 
the model. 

The combined effect of the three fundamental mechanisms of a tur- 
bulent flame, including turbulent eddies, laminar flamelets, and volumetric 
expansion, as discussed in last chapter, elucidate a variety of phenomena 
related to flame stability, effects of confinement and geometrical obstacles, 
flame acceleration and its effect upon flow instability. However, attempts to 
study theoretically these processes and to reproduce the outcome of their 
interaction have been plagued by lack of proper models of turbulent flow, 
difficulties associated with representing chemical reactions in a fluctuating 
field (Mellor & Ferguson, 1980), and numerical instabilities which are intro- 
duced when dynamic effects of combustion are expressed in terms of the 
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expansion across the flame front (Williams, 1974). The crux of the problem 
lies in the effect of the Arrhenius exponential term in the kinetic rate equa- 
tions that describe the chemical reaction. Thus attempts to incorporate 
turbulent fluctuations in the combustion process using statistical decompo- 
sition of thermodynamic variables leads to non-convergent solutions. More- 
over, conventional numerical methods used in these studies are influenced 
by artificial viscosity inhibiting the growth of nature flow instabilities at high 
Reynolds numbers, tending to "laminarize" the flow (McDonald, 1979). On 
the other hand, the use of grids to calculate the flow variables imposes a 
limit on the spatial resolution of the results and may require the added com- 
plicity of adaptive modifications around zones of large gradients with the 
concomitant effects of numerical diffusivity. 

Conventional artificial modeling of turbulence is avoided by the use of 
the random vortex method to evaluate the fluctuating field and the interac- 
tion of turbulent eddies. Our numerical model stresses the unsteady and 
highly fluctuating nature of the flow field associated with turbulent flames, 
as well as the dynamic effects due to exothermicity combustion in a flow sys- 
tem. It was used originally to study the evolution of the vorticity field and 
the development of the flame front behind a step in a relatively short chan- 
nel. Presented in this chapter are solutions we obtained for a long channel 
incorporating a smooth contraction followed by an abrupt expansion. This 
relates to essential geometrical features of the combustion tunnel used for 
experimental studies. In order to treat the appreciably large flow field, we 
had to use a Cray computer with the numerical code properly vectorized for 
this purpose. Besides the vorticity field and flame front kinematics, the 
results include here the detailed average velocity profiles and turbulence 
intensity profiles, in comparison with experimental results of Pitz & Daily 


43 



(1983). 


4.1. Configuration and Computation Parameters 

The overall configuration of the rearward-facing step combustor used 
by Pitz Sc Daily is shown in Figure 14 with all dimensions properly labeled. 
The vorticity field and flame development behind the rearward-facing step, 
without considering any information about inlet section and step geometry, 
was modeled by Ghomem, et al (1982). Their results (see Figure 15) have 
been in agreement with the experimental schlieren cinematographys 
recorded by Ganji Sc Sawyer (1980). 

In order to simulate the reacting flow in greater detail, the computa- 
tional domain is extended to cover the entire step and the upstream inlet 
section. The width and the velocity U- at the inlet of the duct are chosen 
to be the characteristic length and velocity respectively. Since the reat- 
tachment length of the shear layer behind a step is approximately about 6 
to 8 step heights H Q , we here choose the duct length behind the step to be 
10 H q (or 5 H^) to prevent the untractable outflow boundary condition in the 
incomplete recirculation zone. The duct length upstream the step is chosen 
to be 2H^ and the geometry of the step is tried to match the experimental 
shape as closely as possible without introducing too many computational 
complicities. The computational configuration of the rearward-facing step is 
shown in Figure 18. For premixed propane-air reactants with equivalence 
ratio p = 0.57 the density ratio across the flame front is approximately equal 
to 6.0 and the measured laminar flame speed is about 0.02 U Q , where U Q is 
the mean velocity at the step. 
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The computation was performed at two Reynolds numbers Re = 22,000 
and 37,000, based on the inlet flow velocity U^, its width and the kinematic 
viscosity of the reactants v, in order to match the experimental conditions 
of Pitz & Daily. At zero time step, the flow field is initialized by a potential 
flow with a velocity = 1 uniform at the inlet section with a vortex sheet 
layer of thickness 5 = 3 a , where a = (2 t/Re)^ is the standard deviation of 
the set of random displacements mentioned in section 2.2.2, on each wall of 
the channel. The sheet is discretized along each wall into a number of 
finite-length vortex elements, each of h = 0.2 in length. A maximum 
strength r __ = 0.05 for each vortex sheets is chosen to control the error 
in simulating the diffusion process of vorticity; therefore, at each point, 
there will be four sheets generated for each unit wall velocity u w . The dispo- 
sal of vortex sheets and sheets-blobs transformation will then follow the 
algorithm described in our numerical model. 

The Euler flow for this specific geometry was solved by the SOR-vortex 
method. The flow domain was discretized for this purpose into square 
meshes, each of hj = 0.05 in size. The optimal value of relaxation parameter 
was found to be u Q = 1.87. To reach 1 x 10 accuracy of function under 
consideration, the number of iteration for each time step are found to range 
between 20 and 40. 

As vortex sheets are continuously created on the walls, the vorticity 
field of shear layer behind the step and boundary layer along walls start to 
grow; accordingly, the number of vortex blobs increase with time. Vortex 
blobs are disregarded once they move beyond the exit section of the chan- 
nel due to the finite computational domain. The number of vortex elements 
finally oscillate around a constant value and the flow and vorticity field are 
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then considered to reach a stationary state. 

We start our reacting flow calculation when the flow field has been con- 
sidered fully developed. Previously, a point ignition source was provided at 
the corner behind the step to study the development of reacting flow from 
distributed flamelets to flame front (Ghoniem et al, 1982). The number of 
time steps required to establish a continuous flame front is sometimes over 
100. Computational cost of this portion is always expensive, especially when 
the number of vortex elements reaches several thousands. Therefore, 
instead of expressing the initial condition in terms of a point ignition, the 
flame front is introduced as an interface extending along the horizontal 
centerline throughout the flow field in the combustion section to save some 
computational efforts. 

For the sake of convenience, the mesh size used in the flame front cal- 
culation matches that used in non-reacting flow, i.e., h c = hj = 0.05. The 
time step for the reacting flow calculation is then properly reduced to 
satisfy the Courant stability condition. 


4.2. Vorticity Field 

To visualize the flow field all vortex blobs used in the computation are 
plotted as small circles while line segments display vectorial properties of 
their velocities. The developing vorticity fields of Re = 22,000 and Re = 
37,000 are displayed in Figure 17 6c 18 respectively from T = 0.0 to 10.0 for 
every twenty time steps. The vorticity fields exhibit large scale structures 
associated with the generation of finite- size eddies that grow by entraining 
fluid from the inviscid portion of the flow field. The coherence is manifested 
by continuous paring of eddies as they propagate downstream. The 


46 



structure of the mixing layer is located between the free stream of the 
incoming flow and the recirculation zone. For step height equal to one half 
of the channel width, the mixing layer is likely to interact with the boundary 
layer at the upper wall as they both grow downstream. The mixing layer for 
this specific configuration is certainly affected by the presence of the upper 
wall as we can also see in Figure 17 and 18. 


4.3. Flame Development 

Figure 19 depicts the flame front contour and the vorticity field of the 
reacting flow after T = 11.0 for every 0.2 nondimensional time. The flame 
front initialized by a horizontal line across the middle of the combustion 
section, soon follows the flow field and bounds the large scale eddies of the 
mixing layer which provide, through their stirring action, a good contact 
between the products in the recirculation zone and the incoming reactants. 

The vorticity field of reacting flow still clearly exhibits the large scale 
structures which are coherent and grow as they move downstream. This 
suggests that heat release which results in expansion and increase in 
kinematic viscosity of the gas mixture in the mixing layer does not consider- 
ably affect the vortex shedding behind the step. However, the number of 
vortex elements which have a resemblance to eddies in the turbulent flow 
field decreases because of the actions of source blobs. This "laminarization” 
phenomenon in reacting flows (Takagi et al, 1980) results from the increase 
in viscosity of hot products can be modeled without modifying Reynolds 
numbers in the computational parameters. 
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4.4. Haan Velocity Profiles 

Velocity profiles at the step and eight equally-spaced cross sections 
downstream are calculated by averaging the results of 20 consecutive com- 
putational time steps. The vorticity field under consideration is shown in 
Figure 20 and 21. for the cases of Re = 22,000 and 37,000 respectively. 
Since the flow has reached a stationary state and the incoming flow is 
steady, time-averaging is identical to ensemble averaging. While the flow 
field is continuously perturbed by the random samples used in the random 
walk modeling of diffusion, the growth and decay of instability is governed by 
the non-linear interaction due to convection. 

The comparison between our numerical results and experimental data 
of Pitz & Daily for the nonreacting and the reacting cases of both Reynolds 
numbers are displayed in Figure 22 and 23. The results are quite satisfac- 
tory. Since the flow behind the step has a definite separation point and is 
dominated by large eddy structure, which is well reproduced by our method, 
while the fluctuations associated with high Reynolds number flow conditions 
are not encumbered by numerical diffusivity, the agreement between the 
numerical and experimental profiles is especially good in the isothermal 
case. The discrepancies of mean velocity profiles near the step are properly 
due to the fact that the number of vortex blobs in the mixing layer is still 
too small to simulate the large velocity gradient across the layer. This claim 
can be seen clearer in the reacting flow case where the number of vortex 
blobs in the shear layer are reduced by the effect of source blobs, the devia- 
tion of numerical results from experimental data as a result becomes larger. 

Since the reacting flow behind the step is also dominated by the large 
scale structure that causes wrinkling and stretching of the flame front when 
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entraining the reactants into the combustion zone, our numerical model 
also yields quite satisfactory results for mean velocity profiles in the react- 
ing case. The increase of flow velocity in the reaction zone due to thermal 
expansion is properly modeled by volumetric sources along flame front. Our 
phenomenological model of flame propagation has certainly the capability of 
predicting the complicated reacting flow field without solving detailed chem- 
ical kinetics and energy equation. 


4.5. Reattachment Length 

The reattachment length is the length measured from separation point 
of flow field to the reattachment point where the dividing streamline returns 
to the wall. Its value is often expressed in terms of step height and varies as 
a function of some flow parameters such as expansion ratio, Reynolds 
number and boundary layer informations at the step (Eaton and Johnson, 
1901) 

In the non-reacting case, the reattachment length as deduced from the 
mean velocity profiles is about eight step heights for both Reynolds 
numbers, being slightly longer than the data obtained by Pitz k Daily. How- 
ever, when compared with other experimental data collected by Durst and 
Tropea (1982), our results even seem to have closer matches in terms of 
Reynolds number and expansion ratio (Figure 24). 

For reacting flow, the length is decreased to between five and six step 
heights as a consequence of the expansion due to combustion. The recircu- 
lation is also reduced accordingly while the maximum reverse velocity and 
recirculation rate is increased. 
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4.6. Turbulence Intensities 


Turbulence intensities are important parameters in the study of tur- 
bulence. These quantities are difficult to evaluate not only by analytical 
modeling but also by experimental measurements. Even for rearward-facing 
step, there is a substantial variation between experiments in the values of 
the turbulence intensity. The variation is probably caused as much by 
experimental uncertainty as by real difference in the flows. According to 
the review by Eaton & Johnson (1981), the maximum value of streamwise 
turbulence intensity (u'^)^ / U 0 lies from about 0.18 to 0.21 for most exper- 
iments while the peak value of the shear stress -u'v' / U Q is estimated to 
be in the neighborhood of 1.25 x 10 

Figure 25 provides the comparison between the experimental and 
numerical values of streamwise turbulence intensity (u‘ )" / U Q in both 
non-reacting and reacting case. The peak value occurs, as expected from 
the velocity field, at three places : on both walls and at the start of the mix- 
ing layer. Turbulence due to high shear is continuously generated there, 
giving rise to maximum levels, followed by decay downstream. The maximum 
value of (u’ 2 )* / U Q of our numerical results in the shear layer region is 
about 0.18, which is at the lower edge of most experimental data. We also 
have very obviously large intensities on both walls, while the experiment 
data do not show peaks there. The measured intensities also have very nice 
profiles which develop with the shear layer. However, the maximum tur- 
bulence intensity obtained by Pitz & Daily, (u'^)^ / U Q l max = 0.28 .seems to 
overpredict the value by comparing with other measurements. This leads to 
a large discrepancies in the comparison of the maximum value. 
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While numerical results resemble the trend of experimental data, its 
profiles exhibit more irregular shape in the shear layer region and the 
discrepancy becomes particularly apparent downstream. The results also 
show an early decay of the peak value at 2 to 3 step heights upstream the 
reattachment point. All of these deviations might be in part, due to the fol- 
lowing reasons : 

(i) the relatively small size of the sample used in evaluating the numeri- 
cal averages : this is particularly apparent when higher moments of the 
solution are evaluated. Computing them is equivalent to time 
differentiation, emphasizing errors in deviations from the mean. Higher 
moments such as skewness and flatness factor thus were not calculated. 

(ii) the influence of the boundary condition at the exit. This can be 
attributed to two factors : (a) the deletion of vortex blobs at the outflow sec- 
tion to save computational time and (b) the boundary condition d 'f' / d x = 0 
or v = 0 used to evaluate the potential flow velocity. The fluctuating com- 
ponent of the velocity is artificially reduced, as manifested by the fact that 
the comparison is more favorable near the step than towards the end of the 
channel. 

(ill) the less accuracy of the vortex blob method near the boundaries 
where large number of vortex blobs are generated and accumulated. The 
overlapping of vortex blobs cause large fluctuations of streamwise velocity. 

(iv) three-dimensional effect : Experimental results indicate that tur- 
bulence intensity grows around the mixing layer close to the separation 
zone. This is apparently due to three dimensional effects by which energy is 
drawn from the main flow into turbulent fluctuation by vortex stretching, an 
effect beyond the scope of our purely two-dimensional modeL 
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In Figure 28 we present the comparison between the experimental and 
numerical values of transverse turbulence intensity (v‘^)^ / U Q , in both the 
reacting and non-reacting case. Our results have better agreement in this 
direction. The transverse turbulence intensity of numerical calculation is 
shown to have the similar trend and comparable magnitude as (u‘^)^ / U Q . 
while the experimental results show a smaller transverse intensity than the 
streamwise one (Etheridge & Kemp, 1978). 

Since the no-slip boundary condition in our numerical model is satisfied 
by creating new vortices every time step, the random displacements of vor- 
tex elements and the transformation between blobs and sheets in the region 
of high vorticity density decidedly induce considerable fluctuations of 
streamwise velocity. The transverse velocity v, on the other hand, is 
automatically equal to zero when inviscid flow solution is obtained. There- 
fore, the streamwise turbulence intensities obtained by methods based on 
vortex dynamics are always found to have a very large value near walls (Dai 
et al, 1983, Ashurst, 1980) and the transverse turbulence intensities have 
smaller values and converge earlier to zero near walls. 

The discrepancy between experimental turbulence intensities and 
numerical results was also reported by Briggo, et al (1978), who used a Rey- 
nolds stress closure model and a law-of-the-wall boundary condition to com- 
pute the same flow field. Their results either under- or over- predicted the 
streamwise component of turbulence intensity, depending on the rate of 
velocity generation on the step. A recent publication (Walterick, et al, 1984) 
on the numerical calculation of the turbulent flow behind a backward- facing 
step, showed a very good prediction of experimental turbulence intensities. 
They used the k - e model with a special modeling of pressure-velocity corre- 
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lation to match the experimental data. While most of the existing tur- 
bulence model requires artificial modeling of some correlation in the aver- 
aged Navier-Stokes equation, our method can at least qualitatively predict 
turbulent intensities without all these artificial informations. 


4.7. Conclusions 

1. The flow field is dominated by the presence of large scale eddies that 
separate the incoming reactants from the recirculation zone embodying the 
combustion products. The far-field effects of the eddies enhance the 
entrainment of fresh reactants into the mixing zone. 

2. The flame front is stabilized on the contours of the eddies and acts as 
a contact surface between the reactants and the products; as a conse- 
quence, it propagates almost tangentially to the local flow direction. 

3. The average velocity profiles for both the non-reacting and the react- 
ing case are satisfactorily predicted; thus, without modeling the effect of 
turbulent fluctuations, the method provides realistic information about the 
overall dynamic behavior of the flow field. 

4. The prediction of turbulence intensity is less satisfactory because of 
the small size of the sample, and the three-dimensional effects which were 
not included in our calculations. 

5. The length of the recirculation-zone is accurately predicted in both 
the isothermal and the reacting flowand is relatively independent of the 
Reynolds number. 

Our results indicate that the random vortex method is capable of simu- 
lating the large scale, mean behavior of turbulent flow in a two-dimensional 
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field. Turbulent intensity is under-predicted because three-dimensional 
effects are not taken into account. A large ensemble of data as well as a 
precise specification of inlet conditions in terms of vorticity distributions 
are required to accurately evaluate high moments of the velocity field. The 
flame model we adopted was quite adequate to provide information about 
the front geometry as well as the dynamic effect of combustion on the tur- 
bulent flow field. Some of these conclusions were speculated upon by Peters 
and Williams (1980) after Ghoniem, et al published their first paper on the 
model. It is our hope that the new results presented here will promote 
further constructive debate on the subject. 
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Chapter 5 


Unconfined Flame : 

A Turbulent Flame Stabilized by a Circular Flameholder 

Inserting rods in the stream of combustible mixtures is another way of 
stabilizing the flame. The near-wake region behind the rods is capable of 
recirculating combustion products and continuously igniting the fuel mix- 
ture. Fluid mechanics plays an important role in this process. The investi- 
gation of the detailed flow field in the wake region behind the circular flame 
holder can surely give us a better understanding of the mechanisms of flame 
stabilization. 

Two-dimensional turbulent flow behind a circular cylinder has been stu- 
died with the vortex method by Chorin (1973) and Cheer (1980) for only a 
relatively short computational time and Reynolds number up to 10 . Both 
their results have shown very interesting initial flow structures. They also 
predicted accurate drag and lift coefficients with experimental data. 
Another numerical method based on the QUICK method was used by Davis & 
Moore (1982) to study flows behind rectangular cylinders. They obtained a 
series of impressive streakline and isovorticity plots. The computed lift, 
drag and Strouhal number were also found to compare well with experimen- 
tal data for Reynolds number below 1000. 

Recently. Perry, et al (1982) took advantage of the movie made by 
Prandtl to obtain the instantaneous streamline patterns behind a circular 
cylinder for the flow of Reynolds number about 100. A simple model for 
steady vortex shedding was deduced from those streamline patterns. The 
salient features of the flow pattern was also interpreted by using properties 
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of critical points. They conjectured that the vortex shedding process 
behind a circular cylinder is insensitive to Reynolds number. 

The visualization of the flow field in experiments are always accom- 
plished by continuously injecting dyes, gas bubbles or small particles 
upstream the test section. A collection of experimental works of flow visual- 
ization can be found in the album by Van Dyke (1980). Numerically, the 
simulation can be done by releasing a number of passive particles at fixed 
points upstream the body every time step and tracing their positions at con- 
secutive time steps. Streakline plots at any specific time step can be simply 
obtained by displaying the instantaneous locations of all particles which 
have been released. The results of streakline plots obtained earlier by 
Fromm Sc Harlow (1983) and recently by Davis Sc Moore (1982) have attracted 
a considerable amount of attention. 

Although streaklines can be used to give an idea of where the vorticity 
resides in a flow field, they can tell us very little about the surrounding flow 
field and the entrainment process. Instantaneous streamlines, which are 
difficult to be obtained experimentally but are easier to be calculated 
numerically, can serve this need. The particle paths, on the other hand, 
describe the trajectory history of specific particles. It gives us the duration 
of time that particles reside in a given region, thus has application in the 
areas of chemical reaction. 

In this chapter, calculations are extended to a longer computational 
time to study the development of Karman vortex street at high Reynolds 
numbers up to twenty diameters downstream. The instantaneous stream- 
lines in the formation zone are used to study the vortex shedding process 
after the flow is impulsively started. Streaklines and particle pathlines are 
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also computed. 

Flame stabilization is an important process in the practical design of 
combustion chamber. Theories emanating from Spalding (1957) and Zukoski 
fie Marble (1958) in this area are found to correlate data well. Wright fie 
Zukoski (1958) also investigate flame spreading behind flame holders in a 
confined duct and found that flame spreading is independent of free-stream 
velocity, temperature, equivalence ratio and flame speed as long as the flow 
is turbulent and subsonic. 

In an unconflned region, flame stabilized behind a circular cylinder 
always shows a symmetrical wrinkling V-shaped fronts (Petersen fi: Emmons, 
1981, Cheng 8c Ng, 1983 and Yoshida fie Tsuji, 1984). It was found the wrinkles 
have relatively large mean wavelength compared to the integral scales of the 
unburned gas turbulence. Although turbulence intensities of unburned gas 
decrease downstream, the amplitude of wrinkles increases. It has been sug- 
gested (Yoshida fit Tsuji, 1984) that laminar flame instability plays an impor- 
tant role in the mechanism of flame wrinkling. 

Instead of investigating the flame stability and blow-off limit, the results 
presented in this chapter are primarily concerned with the hydrodynamic 
effects on flame shape and burning rates in the near wake region behind a 
circular flame holder. Several different values of burning speed, density 
ratio and Reynolds numbers are used to investigate their effects. 


5.1. Numerical Parameters and Formulations 

For simple geometry like a circular cylinder in our case, the inviscid 
flow component can be easily evaluated by the method of images. The 
potential flow solution of a uniform flow past a cylinder is obtained by the 
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circle theorem (Miline-Thomson, 1970) aa 


W p (z) = 



(5.1) 


where U m is the velocity of the uniform flow at infinity and a is the radius of 
the cylinder. 

In non-reacting flows, for every vortex blob in the flow domain, there is 
a corresponding image vortex of reverse strength at the inverse point a / 
z* to cancel out the normal velocity component at the surface of the 
cylinder (Figure 27). Usually a vortex with equal strength at the center of 
the cylinder is also needed to conserve the vorticity. According to Cheer 
(1980) this vortex should be omitted in order to satisfy the zero total vorti- 
city condition at infinity. The complex velocity induced by all vortex blobs 
can be obtained as 
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In reacting flows, image source blobs have same signs and strengths as 
their corresponding ones along flame front. The complex velocity induced 
by all source blobs can be formulated in a similar way as: 
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(5.3) 
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The total inviscid flow solution thus can be obtained by adding up above 
equations. For the sake of convenience, the diameter d and the uniform 
flow velocity U m are chosen as the characteristic length and velocity for 
nondimensionalization. The time scale is therefore nondimensionalized with 
respect to d / U m . 

To calculate the boundary-layer flow near the surface of the cylinder, 
the surface is divided into two parts : the upper and the lower semi-circle. 
Each part is then discretized into a number of equally spaced segments. 
Since it is easier to compute the flow field in a universal coordinate, each 
semi-circle is transformed into a line segment according to the following 
formula 


x = a 


±tt - tan 1 

x 


(5.4) 


y' - Vx s + y 8 - a (5.5) 

where superscript * denotes the coordinate after transformation ; + and - 
sign in equation (5.4) apply to the upper and the lower semi-circle respec- 
tively. The front stagnation point of the cylinder then becomes the leading 
edge of a flat plate. Vorticities are created and disposed following the algo- 
rithm described in section 2. 1.2.2. After the velocity field in the boundary 
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Layer and the Locations of all vortex sheets are obtained, the coordinate is 
then transformed back to the physical domain in order to check if there are 
any vortex sheets moved from the upper semi-circle to the lower one or vice 
versa. Transformation between vortex sheets and vortex blobs as well as 
locations of newly generated vortex blobs can also be found. 

Since the recirculation zone behind the circular flameholder acts as an 
constant ignition source, the reacting flow calculation is initialized by a line 
ignition source behind the cylinder. A delay of ten computational time steps 
after the flow is impulsively started is used to start the motion of the flame 
so that the flame front can interact with leading vortices. 

The combustion domain behind the cylinder is divided into a number of 
square cells. The size of each cell is chosen to be h c = 0.25 , so that the The 
size of each cell is chosen to be about the same as the distance between 
check points along the surface of the cylinder such that the computational 
time step of flame propagation can match that used for the calculation of 
turbulent flow field. The spatial resolution of flame front is also required to 
be high enough to identify the flame contour without doing any additional 
interpolation. The number of cells in each direction is properly selected to 
ensure that the flame front does not spread out of the computational 
domain. The laminar flame speed normalized with respect to uniform flow 
velocity S u has the value of from 0.02 to 0.04. The density ratio V across 
flame front also varies from 6.0 to 11.0 for any study of their effects on the 

A Q 

flame propagation. Reynolds numbers of 10 and 10 are used in the react- 
ing flow calculation. 
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5.2. Yorticity and Flow Development 

The flow is initially started with a uniform velocity U M = 1. Vorticities 
represented by vortex blobs are gradually generated at walls to satisfy no- 
slip condition and are carried downstream as time goes on. Figure 28. (a)- 
(d) displays the vorticity field in the wake region at Re^ = 10* , where Re^ = 
U M d/ v .from nondimen sional time 1 to 20 for every 10 computational time 
steps. Each point in the figure represents the location of a vortex blob while 
the line segment associate with the point represents the velocity vector of 
the vortex blob as before. These plots of vorticity field not only show the 
vorticity distribution of the turbulent flow, but also give us a streakline like 
picture of the flow structure as we will discuss in the later section. 

The flow first develops symmetrically after impulsively started. It then 
becomes asymmetric when the vortices shed from upper part of the cylinder 
start to interact with the vortices of opposite sign shed from lower part of 
the cylinder. Vortices are then carried downstream by the convection velo- 
city and grow by entraining surrounding inviscid fluid. As the defect velo- 
city in the wake region becomes smaller downstream, the celerity or the 
propagating velocity of each vortex cen-ter gradually increases to the free 
stream velocity. The centers of vortex structures become difficult to iden- 
tify in the figures. However, their positions can be found by subtracting the 
local flow velocity to the free-stream velocity or the predicted celebrity of 
each vortex center (Cantwell Sc Coles, 1983). It can also be noticed from the 
figure that as the flow propagates downstream, the vorticity field spreads 
out wider and wider due to the viscous diffusion (or random walk). The 
number density of vortex blobs inside the large scale structure also 
becomes smaller. This implys that the circulation of large scale stn. ctures 
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decays as the structures move downstream- 


5.3. Vortex- Shedding Process 

For flows at low Reynolds number, the Von Karman vortex street has 
already been recorded by many investigators (Taneda, 1952). It was found 
that the flow pattern still exhibits regular structures up to Reynolds number 
of the order of 10 . However, the detailed vortex shedding process behind 
the cylinder is still not well understood. Recently, Perry et al (1982) made 
an effort to investigate the vortex shedding process for the flow past a 
cylinder or a vertical plate at low Reynolds numbers experimentally. In this 
section, the process will be studied for higher Reynolds numbers by our 
numerical model. 

In order to inquire into the vortex-shedding process, the instantaneous 
streamlines in the vortex formation zone for Reynolds numbers of 10^ and 
10 are calculated. Numerically streamlines can be obtained either by 
integrating the instantaneous velocities 

♦ ( z ) ~ J u dy - - f v dx (5.6) 

or by directly taking the imaginary part of the complex potential 0 (z) = $ 
(z) - i 4' (z). For a uniform non-reacting flow past a circular cylinder, the 
stream function can be formulated as 


* ( * ) = U m y (1 - 
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A continuous streamline is obtained by connecting all the points which have 
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same values of stream function in the flow field. Results in Figure 29 and 30 
show that after the flow is impulsively started, two symmetric separation 
bubbles begin to grow behind the cylinder. As a consequence of the hydro- 
dynamic instability, one of these two bubbles becomes bigger than another 
as shown in Figure 29. a and 30. a. At T = 2. the bigger one starts to detach 
from the surface of the cylinder and moves downstream. The so-called 
"alleyway" forms between two cavities. As the front cavity is gradually car- 
ried downstream, its size becomes smaller and finally opens up as shown in 
Figure 29. b and Figure 30.b. 

Meanwhile, another cavity starts to leave from another part of the sur- 
face and propagates downstream. The sizes of cavities are about the same 
as the diameter of the cylinder for two Reynolds numbers under study. This 
is different from the experimental results of low Reynolds numbers obtained 
by Teneda (1952), for which the size of each cavity can grow to as big as 
three diameters. 

Mostly, cavities will become smaller and finally open up before they 
reach two diameters downstream. Saddle points then cease to exist and 
streamlines become wave-like afterward. It could sometimes happen that 
cavities will entrain more fluids and grow in size again (see Figure 30.b). The 
cavities in this case can remain closed for a longer time. 

By comparing the instantaneous streamlines in Figure 29 and 30 with 
those obtained by others (see, for example, Perry et al, 1984 and Cantwell & 
Coles, 1983), one can find that the vortex shedding process behind a bluff 
body at high Reynolds number has the same qualitative picture as that at 
low Reynolds number except that the vortex formation region reduces in 
size. Separation cavities with two’ different rotating directions shed alterna- 


63 



tively from the back surface of the cylinder. They form centres, saddle 
points and alleyways in the wake region. The centres in the cavities do not 
necessarily match foldings of vortex sheets, while the valleys of the wave- 
like vorticity field always correspond to the saddle points. 


5.4. Streamlines, Pa thlin es and Streaklines 

The unsteady feature of the flow can be examined by the differences 
among streamlines, pathlines and streaklines. A streamline is a space curve 
which is tangent everywhere to a velocity vector, a pathline is the trajectory 
of a fluid particle of fixed identity, while a streakline is the instantaneous 
locus of points that connect the location of all particles which have passed 
through a fixed point in the flow field. They are all identical in a steady flow. 

Instantaneous streamlines at T = 20, as observer stays stationary with 
respect to the cylinder, are shown in Figure 31. a. They exhibit wave-like 
shapes except in the vortex formation zone. Figure 31.b shows the 
corresponding streamlines as observer stays stationary with respect to the 
incoming flow. Centers of vortices do not arrange themselves on two rows 
with a definite spacing ratio as that occurs at lower Reynolds numbers. 
Instead, two neighboring vortices of same rotation tend to get closer, while 
another two neighboring vortices of opposite sign are going apart. We can 
therefore observe two vortices separate another two of opposite sign in an 
alternative way. Two neighboring vortices of same sign also have the ten- 
dency to merge into a single vortex as shown in Figure 31.b. The strong 
interaction of two vortex streets at high Reynolds number was also observed 
by Papailiou and Lykoudis (1974) in their experimental investigation of a 
turbulent wake. 
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Streaklines are calculated by releasing 12 passive particles every com- 
putational time step at 12 equally spaced locations upstream the cylinder. 
Results of the instantaneous locations of all released particles are shown in 
Figure 32 for T = 5.0, 10.0, 15.0 and 20.0 respectively. Streaklines are found 
to be quite chaotic in the wake region, but are continuous and smooth in the 
outer part of the flow. This is consistent with most of the experimental 
results (Van Dyke, 1980). 

By comparing Figure 28 with Figure 32, one finds that the chaotic 
region of streaklines exactly matches the place where vorticities are distri- 
buted. It was also found by Davis & Moore (1982) in their calculation of flows 
behind rectangular cylinders that streaklines exactly match the plots of iso- 
vorticity. This is not surprising because vortices (or eddies) are the funda- 
mental elements of turbulence. 

Figure 33 presents 12 particle pathlines by following their motions from 
time T = 0.0, 5.0, 10.0, 15.0 to T = 20.0 respectively. Consecutive points 
represent their instantaneous positions at the following time step. Dis- 
tances between consecutive points along the same pathline are thus propor- 
tional to their instantaneous velocities. Particles will have larger residence 
time if they have ever been brought into the wake region by the flow field. 


5.5. Flame Development 

Development of turbulent flame behind a circular flameholder are stu- 
died by igniting the flow at the 10th computational time step after the flow is 
started. Results of flame fronts and the corresponding vorticity fields are 
shown in Figure 34. The flame first propagates symmetrically with a slower 
motion near the center line of the cylinder due to the effect of defect 


65 


velocity. At T = 3 the first vortex structure, which just shed from the sur- 
face of the cylinder, starts to stretch the flame front and entrains unburned 
mixtures into the burned region. After some time part of the entrained 
unburned mixtures are surrounded by burned products and form islands 
and pockets inside the flame contours. These unburned islands not only 
increase the burning area, but also widen the flame spreading angle by 
pushing out flame fronts as a result of the volumetric expansion. The burn- 
ing rate is also increased by flame wrinkling and entrainment mechanism of 
turbulent eddies. Computations are performed from non-dimensional time T 
= 0.0 to 15.0. Although the flame front is still in its developing stage, the 
tendency of the flame front to develop into a V-shaped flame can be seen 
from Figure 34. 

Figure 35. a - 35.e provide the comparison of the instantaneous vorticity 
field and velocity vectors between non-reacting and reacting flows from T = 
11.0 to 15.0. Local velocities in the wake region are found to be quite 
fluctuating for both cases. The velocity vectors of reacting flows outside the 
flame boundary indicate the divergence of streamlines. This is in agreement 
with the experimental results reported by Cheng & Ng (1982). The vorticity 
field of the reacting flow has a considerable deterioration of order in the 
vortex sheet and it is also shown to be more diffused than that of the non- 
reacting flow. 

Several other cases of flame development are also studied by varying 
values of normal burning velocity, density ratio and Reynolds number. 
Results of flame fronts at T=15.0 for each case are shown in Figure 36. By 
increasing density ratio and burning speed, flame fronts are found smoother 
and no more unburned pockets exists inside burned region. Figure 37 shows 
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the results of burning rates for each case. Burning rates for all other cases 
are found to be larger initially and then are caught up by the case with 
smaller density ratio and flame speed. This is due to the fact that flame 
area is reduced by the smoothing effect of larger thermal expansion. 
According to Figure 37, Reynolds number has only very little effect on the 
burning rate. But flame front at higher Reynolds number can follow the 
shape of large scale structures better than that at lower Reynolds number 
because of higher turbulence intensities. In our flame model, viscosity 
changes due to combustion are not taken into account. Its value can 
change dramatically due to temperature increase and this is probably one 
of the reasons why symmetrical flame structure is not accurately modeled 
by our numerical method. 


5.6. Conclusions 

1. Periodic vortex street structure exists up to the largest Reynolds 
number under study (lx 10^ ), however, the vortex pairing and diffusion 
process cause some irregularities in vortex spacings. 

2. The size of separation bubble behind the cylinder at high Reynolds 
number is found to be about the same size as the diameter of the cylinder. 
The closed cavity behind the cylinder can sustain only for a short period of 
time. 

3. The size of cavities always reduces downstream and finally opens up. 
However, the cavity may have the chance to entrain more fluids inside itself 
and grows in size. 

4. Vortex shedding process in the formation zone has the same qualita- 
tive picture as long as the Reynolds number is below the critical value. 
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5. Although the flame front is stabilized at the outer contour of eddies, 
the alternating large scale structures are destabilized by the thermal 
expansion of combustion process. The growth of the amplitude of wrinklings 
is associated with the growth of the large coherent structures by the 
diffusion process. 

6. Since the inertial force dominates the flow fields under investigation, 
Reynolds number has little effect on the flame propagation except that the 
flame shape can follow the large scale structures better at higher Reynolds 
number. 

7. Large normal burning velocity and density ratio does not necessarily 
have higher burning rate and spreading angle of flame fronts because they 
tend to smooth out the flame shape and decrease the burning area. On the 
other hand, flame propagation with smaller burning velocity will tend to 
have wrinkler flame front due to stronger interactions of turbulent eddies. 
The unburned pockets inside the flame front , which will push the flame 
more outward when they become burned, can also increase the burning rate 
of the reacting flow. 
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Fig. 1 Configuration of grid points near a curve boundary. 
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Fig. 3 Classification of interface geometries according to SLIC algorithm. 
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Fig. 5 The overall solution procedures for reacting flows. 
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Fig. 6 (a) The coordinate system of a duct with open ends. 

(b) the coordinate system of a duct with one closed end. 
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Fig. 7 


( c ) 

The deformation of flame fronts and velocity vectors in an open 
duct due to three fundmental mechanisms of flame propagation. 
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Fig. B The correponding plots of Fig. 7 in the transformed plane. 
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Fig. 9 The enlargement of the inner portion of Fig. 8. 


77 



Fig. 10 A. phenomenological sketch of flame front stabilized by a circular 
rod holder. 
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Fig. 11 The deformation of flame fronts in a closed duct at the 40th , 120th 
and 200th time step. 
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Fig. 13 The corresponding plots of Fig. 12 in the transformed plane. 
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Fig. 15 The numerical results obtianed by Ghoniem et al (1982) for (a) 
non- reacting flow, and (b) reacting flow behind a rearward-facing 
step. 
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Fig. 17. a 

Fig. 17 The developing vorticity field of the non-reacting flow at Re = 

22 , 000 . 
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Fig. 17. b 
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Fig. 18. a 

Fig. 18 The developing vorticity field of the non-reacting flow at Re = 
37,000. 
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Fig. 18. b 
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Fig. 19 Vorticity fields and flame front contours of the reacting flow at Re = 
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Fig. 20 Vorticity fields at the stationary state of Re = 22,000. 
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Fig. 21 Vorticity fields at the stationary state of Re = 37,000. 
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streamwise velocity profiles for non-reacting and reac 










Fig. 25 Streamwise turbulence intensity profiles for nonreacting and react- 



Fig. 26 Transverse turbulence intensity profiles for nonreacting and react 
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Fig. 20 The developing vorticity field of the wake behind a circular cylinder 
at Re = 10.000. 
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Fig. 28. c 
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Fig. 31 (a) The streamline pattern of Re = 10 4 at T = 20.0 as a uniform flow 
past the cylinder. 

(b) the streamline pattern of Re = 10 4 at T = 20.0 as the cylinder 
moves toward the fluid at a uniform velocity to the left. 
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Fig. 33 Particle pathlines from T = (a) 0, (b) 5, (c) 10. and (d) 15 to T 


respectively. 
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Fig. 34 Flame front contours and vorticity fields of the reacting flow stabil- 
ized by a circular cylinder at Re = 10,000. 
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Fig. 37 Reaction rates of different reacting flows as a function of time. 
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